#include "stdafx.h"
#include "solver.h"

Solver::Solver(double _u0, double _b, double eps, double eps_r, int n, double _h): u0(_u0), b(_b), epsilon(eps), epsilon_right(eps_r), Nmax(n), h(_h)
{

}

double Solver::f(double x, double u)
{
	return u;
}

double Solver::u(double x)
{
	return exp(x);
}

std::list<SearchNode>* Solver::solve()
{
	std::list<SearchNode>* result = new std::list<SearchNode>();

	result->push_back(SearchNode(h,  0.0, u0, 0.0, 0.0, 0.0, u0, u0));

	double x = 0.0;
	double v = u0;
	int n = 0;

	while ((x < (b - epsilon_right)) && (n < Nmax))
	{
		double x_new = x + h;
		double v_new = v + h/2*(f(x,v) + f(x+h,v+h*f(x,v)));

		double v_mid = v + h/4*(f(x,v) + f(x+h/2,v+h/2*f(x,v)));
		double v2 = v_mid + h/4*(f(x + h/2,v_mid) + f(x+h,v_mid+h/2*f(x + h/2,v_mid)));
		double s = (v2 - v_new)*4/3;

		if (abs(s) < epsilon) {
			if (abs(s) > epsilon/8.0)
				result->push_back(SearchNode(h, x_new, v_new, v_mid, v2, s, v_new, u(x_new)));
			else {
				result->push_back(SearchNode(h, x_new, v_new, v_mid, v2, s, v_new, u(x_new)));
				h = h*2;
			}
			x = x_new;
			v = v_new;
		}
		else {
			h = h/2;
		}
		n++;
	}

	return result;
}